# A majority of these codes derive from https://mgimond.github.io/Spatial/interpolation-in-r.html

library(gstat) # Use gstat's idw routine
library(maptools)  # Used for conversion from SPDF to ppp
library(spatstat)  # Used for the dirichlet tessellation function
library(tmap) # Thematic map visualization

# Project sittan_pops_shp

sittan_pops_shp <- SpatialPointsDataFrame(sittan_pops[,c("longitude","latitude")],sittan_pops,proj4string=wgs84) #Transform to spatial vector data frame

region <- spTransform(region,"+proj=eqdc +lon_0=96 +lat_1=16 +lat_2=20 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

sittan_pops_shp <- spTransform(sittan_pops_shp,"+proj=eqdc +lon_0=96 +lat_1=16 +lat_2=20 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

sittan_pops_shp$Total <- sittan_pops$Total

# Replace point boundary extent with that of Bago
sittan_pops_shp@bbox <- region@bbox

# Create an empty grid where n is the Total number of cells
grd              <- as.data.frame(spsample(sittan_pops_shp, "regular", n=50000))
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  # Create SpatialPixel object
fullgrid(grd)    <- TRUE  # Create SpatialGrid object

# Add sittan_pops_shp's projection information to the empty grid
proj4string(grd) <- proj4string(sittan_pops_shp)

# Interpolate the grid cells using a power value of 2 (idp=2.0)
sittan_pops_shp.idw <- gstat::idw(Total ~ 1, sittan_pops_shp, newdata=grd, idp=2.0)

# Convert to raster object then clip to Texas
r       <- raster(sittan_pops_shp.idw)
r.m     <- mask(r, region)

# Plot
tm_shape(r.m) + 
  tm_raster(n=50,palette = "RdBu", auto.palette.mapping = FALSE,
            title="Predicted number of houses") + 
  tm_shape(sittan_pops_shp) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)

# Extract to colonial towns and villages

colonial_towns_shp$pop_1783 <- extract(r.m,colonial_towns_shp,fun=mean) # Population from 1783
sittan_towns_shp$pop_1783 <- extract(r.m,sittan_towns_shp,fun=mean) # Population from 1783

region <- spTransform(region,wgs84)
